#include<stdio.h>
#include<math.h>
//定长矩形法
double rint1(double a, double b, int n, double (*f)(double))
{
	double x, h, s;
	h = (b - a) / n;
	x = a;
	s = 0;
	for (int i = 1; i <= n; i++)
	{
		x = x + h;
		s = s + (*f)(x) * h;
	}
	return s;

}
//变长矩形法
double rint2(double a, double b, double eps, double (*f)(double))
{
    int n = 1;
	double x = a, h, s, s1, p;
	h = (b - a); s = (*f)(x)*h*(*f)(x+h); p = eps + 1.0;
	while (p>=eps)
	{
		s1 = 0.0; x = a;
		h /= 2.0;
        n *= 2;
		for (int i = 1; i <= n; i++)
		{
			x = x + h;
			s1 = s1 + (*f)(x) * h;
		}
		p = fabs(s1 - s);
		s = s1;
	}
	return s;
}
//定长梯形法
double tint1(double a, double b, int n, double (*f)(double))
{

	double x, h, s, fa, fb;

	h = (b - a) / n;
	x = a;
	s = 0;

	for (int i = 0; i < n; i++)
	{
		fa = (*f)(x + i * h); fb = (*f)(x + (i + 1)*h);
		s = s +  h/2*(fa+fb);
	}

	return s;

}
//变长梯形法
double tint2(double a, double b, double eps, double (*f)(double))
{
    double h = (b - a), s = h / 2 * ((*f)(a) + (*f)(b)), p = eps + 1;
    int n = 1;
    while (p >= eps)
    {
        double s1 = 0.0;
        double x = a;
        h /= 2.0;
        n *= 2;
        for (int i = 0; i < n; ++i)
        {
            x = a + i * h;
            s1 += h / 2 * ((*f)(x) + (*f)(x + h));
        }
        p = fabs(s1 - s);
        s = s1;
    }
    return s;
}
//定长辛普森法
double sint1(double a, double b, int n, double (*f)(double))
{

	double x, h, s, fa, fb, fab;

	h = (b - a) / n;
	x = a;
	s = 0;

	for (int i = 0; i < n; i++)
	{
		fa = (*f)(x + i * h); fb = (*f)(x + (i + 1)*h);
		fab = (*f)(x + i * h + h / 2);
		s = s + h / 6 * (fa +4*fab+ fb);
	}
	return s;
}
//变长辛普森法
double sint2(double a, double b, double eps, double (*f)(double))
{
    double h = (b - a), s = 0.0, p = eps + 1;
    int n = 1;
    while (p >= eps)
    {
        double s1 = 0.0;
        double x = a;
        h /= 2.0;
        n *= 2;
        for (int i = 0; i < n; ++i)
        {
            double fa = (*f)(x + i * h);
            double fb = (*f)(x + (i + 1) * h);
            double fab = (*f)(x + i * h + h / 2);
            s1 += h / 6 * (fa + 4 * fab + fb);
        }
        p = fabs(s1 - s);
        s = s1;
    }
    return s;
}
//高斯积分法
 double gaussIntegration(double a, double b, int n, double (*f)(double))
{
    double sum = 0;
    double x, w;
    int i;

    for (i = 0; i < n; i++)
        {
        x = ((b - a) * (i + 0.5) / n) + a; // 计算节点
        w = (b - a) / n; // 计算权重
        sum += w * (*f)(x); // 计算和
    }
    return sum;
}

//差商微分
double CDF(double x, double eps, double (*f)(double))
{

    double h = 0.4;
    double y1=((*f)(x+h)-(*f)(x-h))/(2*h);
    double y2,temp;
    do
    {
        h /= 2;
        y2=((*f)(x+h)-(*f)(x-h))/(2 *h);
        temp = y1;  y1= y2;
    } while(fabs(temp -y2)>eps);

    return y2;
}
//拉格朗日插值求导
double lagrange( double xi,double x[], double y[], int n) {
    double result = 0;
    for (int i = 0; i < n; i++) {
        double term = y[i];
        for (int j = 0; j < n; j++) {
            if (j != i) {
                term *= (xi - x[j]) / (x[i] - x[j]);
            }
        }
        result += term;
    }
    return result;
}

double f1(double x)
{
	return x/(1+x*x);
}
double f2(double x)
{
	return 1+x*x;
}
int main()
{
    double x[3]={1.0,2.0,3.0},y[3]={2.0,5.0,10.0};
	printf("定长矩形法:%e\t变长矩形法:%e\t定长梯形法:%e\t变长梯形法:%e\t定长辛普森法:%e\t变长辛普森法:%e\t高斯积分法:%e\n",rint1(0,3,100,f1),rint2(0,3,1e-6,f1),tint1(0,3,100,f1),tint2(0,3,1e-6,f1),sint1(0,3,100,f1),sint2(0,3,1e-6,f1),gaussIntegration(0,3,100,f1));
	printf("定长矩形法:%e\t变长矩形法:%e\t定长梯形法:%e\t变长梯形法:%e\t定长辛普森法:%e\t变长辛普森法:%e\t高斯积分法:%e\n",rint1(0,2,100,f2),rint2(0,2,1e-6,f2),tint1(0,2,100,f2),tint2(0,2,1e-6,f2),sint1(0,2,100,f2),sint2(0,2,1e-6,f2),gaussIntegration(0,2,100,f2));
	printf("差商微分:%e\t拉格朗日插值求导:%e",CDF(1,1e-6,f2),lagrange(1,x,y,3));

	return 0;

}
